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ABSTRACT 

We simulate the acceleration processes of collisionless particles in a shock structure with mag- 
nctohydrodynamical (MHD) fluctuations. The electromagnetic field is represented as a sum of 
MHD shock solution (B , E ) and torsional Alfven modes spectra (<5B, <5E). We represent fluctu- 
ation modes in logarithmic wavenumber space. Since the electromagnetic fields are represented 
analytically, our simulations can easily cover as large as eight orders of magnitude in resonant 
frequency, and do not suffer from spatial limitations of box size or grid spacing. We determin- 
istically calculate the particle trajectories under the Lorcnz force for time interval of up to ten 
years, with a time step of ~ 0.5 sec. This is sufficient to resolve Larmor frequencies without a 
stochastic treatment. Simulations show that the efficiency of the first order Fermi acceleration 
can be parametrized by the fluctuation amplitude r\ = (SB 2 ) 2 B _1 ■ Convergence of the nu- 
merical results is shown by increasing the number of wave modes in Fourier space while fixing 

V- 

Efficiency of the first order Fermi acceleration has a maximum at r\ ~ 10 1 . The acceleration 
rate depends on the angle between the shock normal and Bo, and is highest when the angle is 
zero. Our method will provide a convenient tool for comparing collisionless turbulence theories 
with, for example, observations of bipolar structure of super nova remnants (SNRs) and shell- like 
synchrotron-radiating structure. 

Subject headings: acceleration of particles — methods: numerical — MHD — turbulence 



Introduction 



Cosmic rays have the spectrum of dN/dE ~ 
10 5 (^/GcV)- 2 - 6 cm- 2 sr- 1 s- 1 GcV _1 up to the so- 
called knee-energy of 10 15 cV. Cosmic ray prop- 
agation theories suggest dN/dE cx E~ 2 energy 
spectra at the cos mic ray acceleration sites (e.g. 



Strong et al.ll2007l ) 



The current description of cosmic ray accelera- 
tion up to knee energy ( 10 15 e V) is the well known 



first-order Fermi acceleration (Axford et al. 1977 



Belli 119781 ). In the first-order Fermi acceleration 
model, magnetic turbulence is an important agent 
for particle acceleration. Turbulence makes the 
particle momenta isotropic, thus portion of the 
particles cross the shock front many times. Expec- 
tation value of the kinetic energy after Nj times of 
shock crossing is E (Nj) = E (1 + h)" J . On the 



other hand, the probability for a particle to sur- 
vive Nj shock crossing can be roughly estimated 
as P (Nj) = (1 — p) Nj . This gives us the power- 
law spectrum of dP/dE = £-(' 1 +p)/'\ 

Ellison et all (Il996t) . iLucek fc Bell (|2000h . and 



Bell fc Lucekl (1200 lb have done simulations to de- 



scribe the self-consistent generation of turbulence, 
with approximations such as gyro-center approxi- 
mation, random walk approximation, or lowering 
the dimension. On the other hand, the recent 
development of the particlc-in cell simulation has 
made it possible to describe the particle acceler- 
ation in_ej^ctron-p_ositron plasma sclf-consistcntly 
(e.g. ISpitkovskvi[2008T) . 

In this letter, we propose an alternative ap- 
proach to the simulation of cosmic ray acccler- 
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ation. We have calculated the motion of par- 
ticles detcrministically, solving the particles' cy- 
clotron motion from Larmor radii of thermal par- 
ticles (~ 10 9 cm) to that of knee energy particles 
(~ 10 17 cm) . According to the theories, we assume 
turbulence spectrum in log k space. This allows us 
to cover a large dynamic-range of space and en- 
ergy, which enables us a direct comparison of the 
accelerated cosmic ray spectra with the observa- 
tions. 

2. Numerical Scheme 

2.1. Representation of Turbulence 

Upstream and Downstream Regions In our 

method, the electromagnetic field and velocity 
field of a continuous region are given by 

B (t, r) = B + Ej B U CX P * ( k i ' r - wjt + <M (1) 
u (t, r) = u + X)j u ij expz (kj • r - uijt + <f>j) (2) 
E(t,r) = -iu(t,r)xB(t,r) (3) 

where the amplitude and the wavenumber of the 
j-th mode are 

Pt 



Table 1: Parameters Used for our Simulations 
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and the initial phase of the j-th mode is (f>j. 

Here Pt is the spectral index that reflects the 
nature of the turbulence and N m is the total num- 
ber of the modes (j G {1, • • • , N m }), Bij is the 
amplitude for each mode, kj is its wavenumber, 
and xii , n 2 are two mutually perpendicular unit 
vectors that are perpendicular to kj . We choose kj 
to be either parallel or antiparallcl to Bq. Equation 

(|6"|) means that ki are logarithmically discrete. 

i 

We use rj = (EjPij 2 ) 2 Bq 1 as the measure 
of the strength of the fluctuation, independent 
of N m . Because increasing N m while keeping 
n = const (1) keeps the magnetic energy in fluc- 
tuation mode, and (2) keeps the expectation value 
of the fluctuation field | (EjBij) | the same, if (f>j 
are independent. We will confirm these properties 
in section [3l 



v Aup 



: 3.0 ■ lCPcms- 1 
= 1.0 ■ 10 7 cms- 1 



fluid speed in shock frame 
Alfvcn speed in fluid frame 



Boup = 1.0-10 5 G unperturbed magnetic field strength 



Jup 

Biup = vBoup 

e 

Amax = 10 17 Cm 
A m in = 10 9 Cm 

T = 0.24kcV 
m = 1.6 • 10" 24 g 
e = 4.8 - 10- 10 esu 



torsional Alfven mode energy 
angle of Bo to shock normal 
maximum wavelength of turbulence 
minimum wavelength of turbulence 
temperature of the particles 
particle mass 
particle charge 



The argument to derive Pt = —1/3 in logfc 
space is summarized below: Variables in log k 
space are marked by tilde. The power law en- 
ergy spectrum is E(k)dk oc k~3 in Kolmogorov 
turbulence case. This energy spectrum is in lin- 
ear bin. In log energy bin the spectral power is 
E(k)d\ogk = kE(k)dk oc k~»\ and since E = 

1/(8tt)B 2 , B(k)d\ogk ex E{kfd\ogk oc k~i. 
Thus, we get Pt = — 1/3 for our discretization of 
the turbulent magnetic field. 

Junction Conditions We assumed strong 
shock junction condition with low plasma f3 limit 



at the shock front: Wd n = J "up, B 



and B 



B 



dn 



B 



|up) 



J, where J is the shock com- 



pression ratio, B\\ and Bi arc components of the 
B normal and tangential to the shock, respec- 
tively. 

2.2. Initial Condition and Equation of Mo- 
tion 

Initial Condition For each set of initial condi- 
tion we introduce electromagnetic fields described 
in section l2~T1 We choose a set of initial turbulence 
phase {4>j}, sign of kj and loj from uniform dis- 
tribution. Then we put 10 5 protons in Boltzmann 
distribution of temperature T at the upstream side 
of the shock. 

We use the va lues in Table [l] based on 
Bamba et all (|2003l )'s observation of SN 1006. 



Evolution We make each turbulence mode 
propagate at Alfvcn velocity of the uniform field 
va — Bo/ y/Airp as in Equation [T] - [3l and up- 
dated the particle with Lorenz force, with 4-th 
order Runge-Kutta method. We choose time dis- 
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Fig. 1. — Particle energy as a function of number of 
shock crossing, after ~ 1 years of time evolution with 

parameters Amax = 10 17 cm , r\ = (HjBij 2 ) 2 Bo _1 = 
10, and 6 = (Bo is parallel to the shock normal). The 
red curves are particle trajectories and inclination of 
the blue line is the prediction of the first order Fermi 
acceleration theory. 
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Fig. 2. — "Convergence" test for energy spectrum. 
Curves shows the particle energy spectrum at ~ 1 year 
time evolution. Each curve corresponds to discretiza- 
tion of the turbulence spectra into log k space with dif- 
ferent Z\log 10 k : number of modes per decade, while 

2\ 



r) = (EjSij 2 ) 5 -Bo -1 = 10 is kept. Other parameters 



are A n 



10 cm and 6 = 0. Particle with energy 



greater than 2GeV are counted. 



at year 1 
at year 3 
at year 1 
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Fig. 3. — Curves show the particle energy spectra for 
Amax = 10 17 cm and rj = 10, 9 = at 1, 3, and 10 years 
of time evolution. Time integrated energy spectrum is 
shown with the bold curve. 



cretization dt for each particle at every timcstcp, 
so that dt < 0.1(1 + r^eB^mT 1 c~ l and dt < 
0.03e |B + E B il m,- 1 ^ 1 always hold. Typical 
time step is 0.5s whereas the Larmor period of 
thermal particle for Bq is ~ 10 2 s. 



dr 
~dt 
dp 
~dt 



P 

e (E 



(e + ^xB) 



(8) 
(9) 



Result 



We have made following examinations for the 
results of our method. First, we have traced the 
particles' energy E as the function of shock cross- 
ing number Nj (Figure [l}. The inclination of 
the curves match the inclinationof the theoretical 
prediction, E(Nj) = {1 + (2/3)(v up - v dn )/c} Nj . 
Secondly, we have traced the spatial location 
where the particles gained their kinetic energy. 
We have found that 94% of the final kinetic en- 
ergy have been earned within 1 final Larmor radius 
away from the shock. This is consistent with the 
first-order Fermi acceleration picture. Thirdly, we 
have studied the validity of our Fourier representa- 
tion in log k space. We have kept the physical pa- 
rameters and increased the number of modes per 
decade Alog 10 k = (N modc / log 10 )(fc max /fc min ); we 
see that the spectra converge, and converge to the 
theoretical power-law spectrum (Figure [2]). This 
justifies our use of log k space discretization. 

We have done a large number of simulations 
while varying the background conditions, A max 
from 10 13 cm to 10 17 cm, 77, the ratio of magnetic 
energy in fluctuation mode to that in background 
field from 0.3 to 300, 6, the angle between the 
background field and the shock normal from to 
tt/2. Figure [3] shows the time evolution of the 
energy spectrum for 10 years. The high-energy 
end of the spectrum gradually grows, and reaches 
2.5 x 10 13 eV by 10 years. 

In our simulation all the particles start its mo- 
tion in the given time. Since we don't include 
the back-reaction from the particles to the electro- 
magnetic field in our simulations, time-integral of 
energy spectra at each time-slice gives the steady 
state energy spectra. This steady state spectrum 
is also shown in Figure [3] with thick curve. The 
nonthermal component has E(k) oc k~ 16 power- 



law spectrum that meets the observational re- 
quirement mentioned in section [1] We can also 
estimate the "injection rate" to be the proportion 
of particles that have more than 2GeV of energy 
after 1 years. For T = 0.24keV, 24keV, 2.4MeV, 
and 0.24GeV, the injection rate was < 0.001, 
1.9 x 10~ 2 , 8.4 x 10 -2 , and 0.378, respectively.The 
other parameters are r\ = 30, 9 = and A max = 
10 17 cm. 

In Figure 2] we show for all the parameter range 
the ratio of the particle numbers that were accel- 
erated to have energy greater than 2GeV. We see 
that the acceleration is most efficient at polar re- 
gion (6 ~ 0) when 77 > 1. We can understand 
this dependence of the spectra with background 
fluid parameters as follows; particles are trapped 
in Larmor motion and tend to move in direction 
of B . Thus particles more easily cross the shock- 
front when B is parallel to shock normal. If the 
turbulence amplitude is much weaker, fewer par- 
ticles get reflected by pitch angle scattering, and 
Fermi acceleration is suppressed. The injection is 
more efficient for smaller A max , because more en- 
ergy is distributed to modes with smallest wave- 
lengths which are resonant with the thermal par- 
ticles. 

If a spherical shock emerges in a uniform mean 
magnetic field, there are two polar region where 
the mean magnetic field is parallel to the shock 
normal, and the equatorial region has the mean 
magnetic field perpendicular to the shock normal. 
Thus, we expect the Fermi acceleration process to 
be only active in the pair of polar region. This 
might explain the bipolar structure we see at SN 
1006. 

We have also checked the acceleration rate 
in three-dimcnsional(isotropic), rather than one- 
dimcnsional(anisotropic) distribution of k^. Wc 
have found that less significant dependence of in- 
jection rate on 6 with larger r\. We can inter- 
pret this as follows; if the turbulence spectrum is 
isotropic and the maximum turbulence wavelength 
is large, turbulence modes with largest wavelength 
and strongest amplitude play the role of local B . 
Thus we observe almost isotropic Fermi accelera- 
tion. This might explain the many SNRs with no 
typical orientation. 



Particles accelerated over 2GeV 




Fig. 4. — Parameter dependence of particle accel- 
eration efficiency . A max is the longest wavelength 

of the turbulence modes, r\ — (T,jBij 2 -Bo -2 ) 2 £ 
{0.3,1,3,30,300} is the ratio of turbulent magnetic 
field to unperturbed magnetic field, 9 is the angle be- 
tween shock normal and Bq. 



4. Discussion 

Some might question the vali dity of 77 value 
much greater than unity. However. lUchivama et al 
(|2007l ) observed extremely fast varying X-ray im- 
ages at SNR RX J1713. 7-3946. Their observation 
may indicate that magnetic field is locally en- 
hanced up to lmG in ~ 1 year in SNR, which cor- 
responds to 77 = 100 case in our model. Our sim- 
ulations suggest that such fast-variating spots in 
SNRs might be the sites of galactic (E < 10 15 eV) 
cosmic ray acceleration. 

Although we have ignored many of the Fourier 
modes by adopting log k space discretization of the 
turbulence spectrum, the validity of the approxi- 
mation can be argued in many ways. Most impor- 
tantly we have confirmed that our measure in log k 
space lead to convergence. Figure [2] shows energy 
spectra for 77 = const, with increasing Z\log 10 k. 
Turbulent cascade, from which the very turbu- 
lence arises, is by nature a logarithmic process: 
a mode of a certain wavelength couples with the 
mode of half the wavelength by nonlinear term 
of Euler equation. First order Fermi acceleration 
also is a logarithmic process: particles gain en- 
ergy as an exponential function of shock crossing 
number E (Nj) = E (l + h) Nj . All these rea- 
son combined, waves in logarithmically discretized 
wavenumber space act as a sufficient ladder to 
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carry up cosmic ray particles. 
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